## VAA effects meta-analysis ----------------------
## sensitivity checks -----------------------------

# load packages, data, and functions
source("packages.R")
load("vaa_df_prep.RData")
source("functions.R")

############################################################################
## Figure C1. Trim-and-Fill: Publication bias ------------------------------
############################################################################

# turnout
turnout_raneff <- rma(yi = log_odds, 
                      sei = log_odds_se, 
                      data = vaa_turnout,
                      slab = study_label,
                      method = "REML")
summary(turnout_raneff)

# vote choice
vote_raneff <- rma(yi = log_odds, 
                   sei = log_odds_se, 
                   data = vaa_vote,
                   slab = study_label,
                   method = "REML")
summary(vote_raneff)

# Figure C1. Trim-and-fill funnel plots random-effects models
pdf(file = "trimfill-turnout.pdf", height = 5, width = 7, family = "Times")
par(oma = c(0,0,0,0)+.1)
par(mar = c(4,4,1,1)+.1)
turnout_raneff_trimfill <- trimfill(turnout_raneff)
funnel(turnout_raneff_trimfill, main = "", legend = TRUE)
dev.off()

pdf(file = "trimfill-vote.pdf", height = 5, width = 7, family = "Times")
par(oma = c(0,0,0,0)+.1)
par(mar = c(4,4,1,1)+.1)
vote_raneff_trimfill <- trimfill(vote_raneff)
funnel(vote_raneff_trimfill, main = "", legend = TRUE)
dev.off()



############################################################################
## Figure B8. Sensitivity to model specification by outcome ----------------
############################################################################

## CCRE models -------------

# turnout
turnout.ccrem <- rma.mv(yi = log_odds, 
                        V = log_odds_se^2, 
                        random = list(~ 1 | study_label, ~ 1 | election),
                        data = vaa_turnout,
                        slab = study_label,
                        method = "REML",
                        tdist = TRUE)

# vote
vote.ccrem <- rma.mv(yi = log_odds, 
                     V = log_odds_se^2, 
                     random = list(~ 1 | study_label, ~ 1 | election),
                     data = vaa_vote,
                     slab = study_label,
                     method = "REML")


## FE models -------------

# turnout
turnout.fixeff <- rma(yi = log_odds, 
                      sei = log_odds_se, 
                      data = vaa_turnout,
                      slab = study_label,
                      method = "FE")
summary(turnout.fixeff)

# vote choice
vote.fixeff <- rma(yi = log_odds, 
                   sei = log_odds_se, 
                   data = vaa_vote,
                   slab = study_label,
                   method = "FE")
summary(vote.fixeff)

# political knowledge
knowledge.fixeff <- rma(yi = pcor, 
                        sei = pcor_se, 
                        data = vaa_knowledge,
                        slab = study_label,
                        method = "FE")
summary(knowledge.fixeff)


## RE models -------------

# turnout
turnout.raneff <- rma(yi = log_odds, 
                      sei = log_odds_se, 
                      data = vaa_turnout,
                      slab = study_label,
                      method = "REML")
summary(turnout.raneff)

# vote choice
vote.raneff <- rma(yi = log_odds, 
                   sei = log_odds_se, 
                   data = vaa_vote,
                   slab = study_label,
                   method = "REML")
summary(vote.raneff)

# political knowledge
knowledge.raneff <- rma(yi = pcor, 
                        sei = pcor_se, 
                        data = vaa_knowledge,
                        slab = study_label,
                        method = "REML")
summary(knowledge.raneff)



## RE models, no weighting -------------

# turnout
turnout.raneff.noweights <- rma(yi = log_odds, 
                                sei = log_odds_se, 
                                data = vaa_turnout,
                                slab = study_label,
                                method = "REML",
                                weighted = FALSE)

# vote choice
vote.raneff.noweights <- rma(yi = log_odds, 
                             sei = log_odds_se, 
                             data = vaa_vote,
                             slab = study_label,
                             method = "REML",
                             weighted = FALSE)

# political knowledge
knowledge.raneff.noweights <- rma(yi = pcor, 
                                  sei = pcor_se, 
                                  data = vaa_knowledge,
                                  slab = study_label,
                                  method = "REML",
                                  weighted = FALSE)


## RE models, giving every study equal weight -------------

# turnout
n_studies <- length(unique(vaa_turnout$study_label))
n_effects_per_study <- as.data.frame(table(vaa_turnout$study_label), stringsAsFactors = FALSE)
names(n_effects_per_study) <- c("study_label", "n_effects")
vaa_turnout <- merge(vaa_turnout, n_effects_per_study, by = "study_label", all.x = TRUE)
vaa_turnout$study_weight <- 1/n_studies * 1/vaa_turnout$n_effects

turnout.raneff.studyweights <- rma(yi = log_odds, 
                                   sei = log_odds_se, 
                                   data = vaa_turnout,
                                   slab = study_label,
                                   method = "REML",
                                   weights = study_weight)

# vote choice
n_studies <- length(unique(vaa_vote$study_label))
n_effects_per_study <- as.data.frame(table(vaa_vote$study_label), stringsAsFactors = FALSE)
names(n_effects_per_study) <- c("study_label", "n_effects")
vaa_vote <- merge(vaa_vote, n_effects_per_study, by = "study_label", all.x = TRUE)
vaa_vote$study_weight <- 1/n_studies * 1/vaa_vote$n_effects

vote.raneff.studyweights <- rma(yi = log_odds, 
                                sei = log_odds_se, 
                                data = vaa_vote,
                                slab = study_label,
                                method = "REML",
                                weights = study_weight)

# political knowledge
n_studies <- length(unique(vaa_knowledge$study_label))
n_effects_per_study <- as.data.frame(table(vaa_knowledge$study_label), stringsAsFactors = FALSE)
names(n_effects_per_study) <- c("study_label", "n_effects")
vaa_knowledge <- merge(vaa_knowledge, n_effects_per_study, by = "study_label", all.x = TRUE)
vaa_knowledge$study_weight <- 1/n_studies * 1/vaa_knowledge$n_effects

knowledge.raneff.studyweights <- rma(yi = pcor, 
                                     sei = pcor_se, 
                                     data = vaa_knowledge,
                                     slab = study_label,
                                     method = "REML",
                                     weights = study_weight)




## CCREM, without data duplicates -----------------------------

vaa_turnout_dup <- vaa_turnout %>% filter(!(election %in% c("Finland 2003", "Germany 2009", "Netherlands 2006", "Switzerland 2007", "Switzerland 2011") & 
                                              bibkey %in% c("Garzia2017", "Marschall2012", "Ruusuvirta2009", "Gemenis2014", "Germann2019", "Germann2019")))

vaa_vote_dup <- vaa_vote %>% filter(!(election %in% c("Netherlands 2006", "Switzerland 2011") & bibkey == "Andreadis2014" |
                                        election == "Netherlands 2010" & bibkey == "Klein-Kranenburg2015"))


# turnout without duplicates
turnout.ccrem.nodups <- rma.mv(yi = log_odds, 
                               V = log_odds_se^2, # square standard errors to feed sampling variances
                               random = list(~ 1 | study_label, ~ 1 | election),
                               data = vaa_turnout_dup,
                               slab = study_label,
                               method = "REML",
                               tdist = TRUE)
summary(turnout.ccrem.nodups)
compute_I2(data = vaa_turnout_dup, model = turnout.ccrem.nodups, vi = "log_odds_se")

# vote choice without duplicates
vote.ccrem.nodups <- rma.mv(yi = log_odds, 
                            V = log_odds_se^2, # square standard errors to feed sampling variances
                            random = list(~ 1 | study_label, ~ 1 | election),
                            data = vaa_vote_dup,
                            slab = study_label,
                            method = "REML")
summary(vote.ccrem.nodups)
compute_I2(data = vaa_vote_dup, model = vote.ccrem.nodups, vi = "log_odds_se")


# collect results

dat.estimates <- data.frame(
  model = c("CCREM","CCREM, no dups", "FE", "RE", "RE, equal weights", "RE, study weights"),
  turnout.b = c(turnout.ccrem$b, turnout.ccrem.nodups$b, turnout.fixeff$b, turnout.raneff$b, turnout.raneff.noweights$b, turnout.raneff.studyweights$b),
  turnout.se = c(turnout.ccrem$se, turnout.ccrem.nodups$se, turnout.fixeff$se, turnout.raneff$se, turnout.raneff.noweights$se, turnout.raneff.studyweights$se),
  vote.b = c(vote.ccrem$b, vote.ccrem.nodups$b, vote.fixeff$b, vote.raneff$b, vote.raneff.noweights$b, vote.raneff.studyweights$b),
  vote.se = c(vote.ccrem$se, vote.ccrem.nodups$se, vote.fixeff$se, vote.raneff$se, vote.raneff.noweights$se, vote.raneff.studyweights$se),
  know.b = c(NA, NA, knowledge.fixeff$b, knowledge.raneff$b, knowledge.raneff.noweights$b, knowledge.raneff.studyweights$b),
  know.se = c(NA, NA, knowledge.fixeff$se, knowledge.raneff$se, knowledge.raneff.noweights$se, knowledge.raneff.studyweights$se),
  stringsAsFactors = FALSE)

dat.estimates$turnout.95lo <- dat.estimates$turnout.b - 1.96 * dat.estimates$turnout.se
dat.estimates$turnout.95hi <- dat.estimates$turnout.b + 1.96 * dat.estimates$turnout.se
dat.estimates$vote.95lo <- dat.estimates$vote.b - 1.96 * dat.estimates$vote.se
dat.estimates$vote.95hi <- dat.estimates$vote.b + 1.96 * dat.estimates$vote.se
dat.estimates$know.95lo <- dat.estimates$know.b - 1.96 * dat.estimates$know.se
dat.estimates$know.95hi <- dat.estimates$know.b + 1.96 * dat.estimates$know.se


# build summary coefficient plot
dat.estimates.rev <- map_df(dat.estimates, rev) # reverse

## create plot
pdf(file="model-sensitivity.pdf", height=3.5, width=8, family="Times")
par(oma=c(1,0,.5,0))
par(mar=c(3,8,2,1))
layout(matrix(1:3, 1, byrow = TRUE), widths=c(1.4, 1, 1), heights=c(2, 2, 2))
colors <- c("black", "darkgrey")

## turnout panel
plot(0, 0, xlim = c(0, 4), ylim = c(.7, 6.3), pch = 20, col = "white", xlab = "", ylab = "", yaxt = "n", xaxt = "n")
axis(2, at = 0:3, labels = F, tick = F)
text(y = 1:6, par("usr")[1], labels = dat.estimates.rev$model, srt = 0, pos = 2, xpd = TRUE, cex = 1)
axis(1, at = seq(0, 4, 1), labels = seq(0, 4, 1), tick = T)
mtext("Turnout", side = 3, line = .7, outer = FALSE)
abline(v = 1, lty = 2)
abline(h = c(.5, 1.5, 2.5, 3.5, 4.5, 5.5), lty = 1, col = "darkgrey")
for(i in 1:nrow(dat.estimates.rev)){
  arrows(x0 = exp(dat.estimates.rev$turnout.95lo)[i], y0 = i, x1 = exp(dat.estimates.rev$turnout.95hi)[i], y1 = i, length = .02, angle = 90, code = 3, lwd = 1, col = colors[1])
  points(exp(dat.estimates.rev$turnout.b[i]), i)
}

## vote choice panel
par(mar=c(3,0,2,1))
plot(0, 0, xlim = c(0, 4), ylim = c(.7, 6.3), pch = 20, col = "white", xlab = "", ylab = "", yaxt = "n", xaxt = "n")
axis(2, at = 0:3, labels = F, tick = F)
axis(1, at = seq(0, 4, 1), labels = seq(0, 4, 1), tick = T)
mtext("Vote choice", side = 3, line = .7, outer = FALSE)
abline(v = 1, lty = 2)
abline(h = c(.5, 1.5, 2.5, 3.5, 4.5, 5.5), lty = 1, col = "darkgrey")
for(i in 1:nrow(dat.estimates.rev)){
  arrows(x0 = exp(dat.estimates.rev$vote.95lo)[i], y0 = i, x1 = exp(dat.estimates.rev$vote.95hi)[i], y1 = i, length = .02, angle = 90, code = 3, lwd = 1, col = colors[1])
  points(exp(dat.estimates.rev$vote.b[i]), i)
}

## knowledge panel
par(mar=c(3,0,2,1))
plot(0, 0, xlim = c(-.05, .25), ylim = c(.7, 6.3), pch = 20, col = "white", xlab = "", ylab = "", yaxt = "n", xaxt = "n")
axis(2, at = 0:3, labels = F, tick = F)
axis(1, at = seq(-.05, .25, .05), labels = seq(-.05, .25, .05), tick = T)
mtext("Issue knowledge", side = 3, line = .7, outer = FALSE)
abline(v = 0, lty = 2)
abline(h = c(.5, 1.5, 2.5, 3.5, 4.5, 5.5), lty = 1, col = "darkgrey")
for(i in 1:nrow(dat.estimates.rev)){
  arrows(x0 = dat.estimates.rev$know.95lo[i], y0 = i, x1 = dat.estimates.rev$know.95hi[i], y1 = i, length = .02, angle = 90, code = 3, lwd = 1, col = colors[1])
  points(dat.estimates.rev$know.b[i], i)
}

# joint x axis label
mtext("Overall effect estimate", side = 1, line = -.5, outer = TRUE, at = .58, cex = 1)
dev.off()


############################################################################
## Figure B9- GOSH plots (graphical display of heterogeneity) --------------
############################################################################

## GOSH plots -------------

# fit FE models to all possible subsets (timely, alternatively load models)
#load(gosh-sims.RData)
turnout.fixeff.gosh <- gosh(turnout.fixeff)
vote.fixeff.gosh <- gosh(vote.fixeff)

# export simulation data
#save(turnout.fixeff.gosh, vote.fixeff.gosh, file = "gosh-sims.RData")

# create GOSH plot
png(file="gosh-turnout.png",
    width = 1000, height = 500)
plot(turnout.fixeff.gosh, breaks = 100, hh = .1, cex.axis = 2, cex.lab = 2, labels = c("Overall estimate (log odds ratio)", expression(I^2)))
dev.off()

png(file="gosh-vote.png",
    width= 1000, height= 500)
plot(vote.fixeff.gosh, breaks = 100, hh = .1, cex.axis = 2, cex.lab = 2, labels = c("Overall estimate (log odds ratio)", expression(I^2)))
dev.off()

############################################################################
## Forest plots for CCREM models without data duplicates -------------------
############################################################################

## forest plots for CCREM without duplicates -------------

pdf(file = "forest-turnout-dup.pdf", height = 5, width = 8, family = "URWTimes")
par(oma = c(1,0,0.5,0))
par(mar = c(1.5,3,0.5,.5))
header_line <- nrow(vaa_turnout_dup) + 2
forest.rma.mod(turnout.ccrem.nodups,
               slab = vaa_turnout_dup$study_label,
               xlim = c(-12, 6), 
               at = log(c(0.25, 1, 5, 12)), 
               atransf = exp, 
               showweights = TRUE,
               ilab = cbind(vaa_turnout_dup$country, vaa_turnout_dup$election_year, vaa_turnout_dup$sample_size, vaa_turnout_dup$vaa_usage_sample), 
               ilab.xpos = c(-8, -5, -3.5, -2), 
               ilab.pos = c(4, 2, 2, 2),
               cex = 0.75,
               xlab = "",
               mlab = "Overall estimate (CCREM with study and election REs)")
op <- par(cex = 0.75, font = 2) 
text(-12, header_line-.2, "Author(s) and Year", pos = 4) 
text(c(-8, -5, -3.5, -2, 1.5, 3.8, 6),  header_line-.2, c("Country", "Year", "Size", "VAA users", "Observed outcome", "Weight", "OR [95% CI]"), pos = c(4, 2, 2, 2, 2, 2, 2))
text(c(-6.5, -3.25, 2.25), header_line + 2.5, c("Election", "Sample", "Estimates"), pos = c(1, 1))
arrows(x0 = -7.75, x1 = -5.25, y0 = header_line + .75, length = 0)
arrows(x0 = -4.25, x1 = -2.25, y0 = header_line + .75, length = 0)
arrows(x0 = -1, x1 = 5.75, y0 = header_line + .75, length = 0)
par(op)
dev.off()

pdf(file = "forest-vote-dup.pdf", height = 4.2, width = 8, family = "URWTimes")
par(oma = c(1,0,.5,0))
par(mar = c(1.5,3,0.5,.5))
header_line <- nrow(vaa_vote_dup) + 2
forest.rma.mod(vote.ccrem.nodups,
               slab = vaa_vote_dup$study_label,
               xlim = c(-12, 6), 
               at = log(c(0.25, 1, 5, 12)), 
               atransf = exp, 
               showweights = TRUE,
               ilab = cbind(vaa_vote_dup$country, vaa_vote_dup$election_year, vaa_vote_dup$sample_size, vaa_vote_dup$vaa_usage_sample), 
               ilab.xpos = c(-8, -5, -3.5, -2), 
               ilab.pos = c(4, 2, 2, 2),
               cex = 0.75,
               xlab = "",
               mlab = "Overall estimate (CCREM with study and election REs)")
op <- par(cex = 0.75, font = 2) 
text(-12, header_line-.2, "Author(s) and Year", pos = 4) 
text(c(-8, -5, -3.5, -2, 1.5, 3.8, 6),  header_line-.2, c("Country", "Year", "Size", "VAA users", "Observed outcome", "Weight", "OR [95% CI]"), pos = c(4, 2, 2, 2, 2, 2, 2))
text(c(-6.5, -3.25, 2.25), header_line + 2.5, c("Election", "Sample", "Estimates"), pos = c(1, 1))
arrows(x0 = -7.75, x1 = -5.25, y0 = header_line + .75, length = 0)
arrows(x0 = -4.25, x1 = -2.25, y0 = header_line + .75, length = 0)
arrows(x0 = -1, x1 = 5.75, y0 = header_line + .75, length = 0)
par(op)
dev.off()
